AEDCTR-87-36 


MAR  2  1 1988 


User’s  Manual  for  EULER3DS 


R.  G.  Hindman 
Iowa  State  University 
Ames,  lA  SOOlO 

PROPERTY  OF  U.S.  AIR  FORCE 
AEDC  TECHNICAL  LIBRARY 

January  1988 


Final  Report  for  Period  October  1 986  -  September  1 987 

technical  reports 

file  copy 


Approved  for  public  release;  distn button  unlimited. 


ARNOLD  ENGINEERING  DEVELOPMENT  CENTER 
ARNOLD  AIR  FORCE  RASE,  TENNESSEE 
AIR  FORCE  SYSTEMS  COMMAND 
UNITED  STATES  AIR  FORCE 


NOTICES 


When  U.  S.  Government  drawings,  specifications,  or  other  data  are  used  for  any  purpose  other  than 
a  definitely  related  Government  procurement  operation,  the  Government  thereby  incurs  no  responsibility 
nor  any  obligation  whatsoever »  and  (he  fact  that  the  Government  may  have  formulated,  furnished,  or 
in  any  way  supplied  the  said  drawings,  specifications,  or  other  data,  is  not  to  be  regarded  by  implication 
or  otherwise,  or  in  any  manner  licensing  the  holder  or  any  other  person  or  corporation,  or  conveying 
any  rights  or  permission  to  manufacture,  use,  or  sell  any  patented  invention  that  may  in  any  way  be 
related  thereto. 

Qualified  users  may  obtain  copies  of  this  report  from  the  Defense  Technical  Information  Center. 


References  lo  named  commercial  producLs  in  this  report  are  not  to  be  considered  in  any  sense  as  an 
endorsement  of  the  product  by  the  United  States  Air  Force  or  the  Government. 


This  report  has  been  reviewed  by  the  Office  of  Public  Affairs  (PA>  and  is  releasable  to  the  National 
Technical  information  Service  (NTIS).  At  NTIS,  it  will  be  available  to  the  general  public,  including  foreign 
nations. 


APPROVAL  STATEMENT 


This  report  has  been  reviewed  and  approved. 


l^-xddL  ^ 


KEITH  L.  KUSHMAN 
Chief,  Facility  Technology  Division 
Directorate  of  Technology 
Deputy  for  Operations 


Approved  for  publication: 
FOR  THE  COMMANDER 


MARION  L.  LASTER 
Technical  Director 
Directorate  of  Technology 
Deputy  for  Operations 


ECURITY  CLASSIFICATION  OF  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


lb.  RESTRICTIVE  MARKINGS 


1a.  REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


2b.  DECLASSIFICATION /DOWNGRADING  SCHEDULE 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBER($) 

AEDC-TR-87-36 


3.  DISTRIBUTION /AVAILABILITY  OF  REPORT 

Approved  for  public  release;  distribution 
is  unlimited. 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


6a.  NAME  OF  PERFORMING  ORGANIZATION  6b.  OFFICE  SYMBOL  7a.  NAME  OF  MONITORING  ORGANIZATION 

R.G,  Hindman,  Consultant  (if applicable) 


6c.  ADDRESS  {City,  State,  and  ZIP  Code)  7b,  ADDRESS  (C/ty,  State,  and  ZIP  Code) 

Iowa  State  University 
Ames,  lA  50010 


8a.  NAME  OF  FUNDING /SPONSORING  8b,  OFFICE  SYMBOL  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

ORGANIZATION  /^j^nold  Engineering  Of  applicable)  rrciQ/i  r 


Development  Center 


8c.  ADDRESS  (C/fy,  State,  and  ZIP  Code) 

Air  Force  Systems  Command 

Arnold  Air  Force  Base,  TN  37389-5000 


1 1 .  TITLE  (Include  Security  Classification) 

User's  Manual  for  EULER3DS 


12.  PERSONAL  AUTHOR(S) 

Hindman,  R.G.,  Consultant,  Iowa  State  University 


CFSI84-5 


10.  SOURCE  OF  FUNDING  NUMBERS 


O. 


PROJECT 

NO. 

DB97VW 


WORK  UNIT 
ACCESSION  NO. 


13a.  TYPE  OF  REPORT 

Final 


16.  SUPPLEMENTARY  NOTATION 


13b.  TIME  COVERED  114.  DATE  OF  REPORT  (Year,  Month,  Day)  115.  PAGE  COUNT 

FROM  10/86  TO  9/87  I  January  1988  I  34 


Available  in  Defense  Technical  Information  Center  (DTIC). 


FIELD 

GROUP 

20 

04 

COSATI  CODES 


SUB-GROUP 


18.  SUBJECT  TERMS  {Continue  on  reverse  if  necessary  and  identify  by  block  number) 

three-dimensional  space  marching  code 
steady  Euler  equations 
supersonic  external  flows 


1 9.  ABSTRACT  {Continue  on  reverse  if  necessary  and  identify  by  block  number) 

Theoretical  background  is  presented  for  a  three-dimensional  space  marching  code  used 
to  solve  the  steady  Euler  equations  for  supersonic  external  flows.  The  ALPHA  scheme  is 
introduced  as  a  novel  second-order  integration  scheme.  The  code  is  coupled  to  the  QUICK- 
Geometry  modeling  system  and  incorporated  either  a  perfect  gas  or  equilibrium  real  gas  mod 
el.  A  discussion  of  numerical  aspects  required  for  a  user  to  set  up  the  grid  and  the  ini¬ 
tial  solution  is  presented. 


20.  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT  21.  ABSTRACT  SECURITY  CLASSIFICATION 

□  unclassified/unlimited  [g  SAME  AS  RPT.  □  pTic  USERS  Unclassified  _ _ 

22a,  NAME  OF  RESPONSIBLE  INDIVIDUAL  22b7TELEPKONE7i^^  Area  Coderi^.  OFFIcTsyMBO^^ 

C.L.  Garner  (615)  454-7813  DOCS 


DD  Form  1473,  JUN  86  Previous  editions  are  obsolete.  SECURITY  CLASSIFICATION  OF  THIS  PAGE 

UNCLASSIFIED 


AEDC-TR-87-36 


PREFACE 

The  work  reported  herein  was  conducted  by  the  author  as  principal  investigator  on  subcontract  CFSI84-5 
between  Iowa  State  University  and  Calspan  Corporation/ AEDC  Division,  operating  contractor  for  the 
Aerospace  Flight  Dynamics  testing  effort  at  Arnold  Engineering  Development  Center,  Air  Force  Systems  Com¬ 
mand,  Arnold  Air  Force  Base,  Tennessee,  and  on  subcontracts  CFSI85-7  and  CFSI86-02  between  The  University 
of  Texas  at  Arlington  and  Calspan  Corporation/ AEDC  Division. 
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A  User’s  Manual  For 
EULER3DS 

by 

Richard  G.  Hindman 

Associate  Professor  of  Aerospace  Engineering 
Iowa  State  University 


1.0  INTRODUCTION 


EULER3DS  is  a  three  dimensional  space  marching  code  for  computing  steady  state  solutions  to 
the  Euler  equations  applied  to  external  flows.  The  code  solves  the  integral  equations  of  motion  in  a 
finite  volume  sense  and  is  thus  totally  conservative.  The  algorithm  used  is  called  the  ALPHA’-scheme 
and  was  developed  by  the  author  for  this  code.  This  scheme  is  explicit  with  either  first  or  second 
order  accuracy  and  is  constructed  by  splitting  the  fluxes  on  cell  faces  to  account  for  upwind  signal 
propagation.  The  code  solves  the  governing  equations  using  either  a  perfect  gas  model  or  an  equilib¬ 
rium  air  gas  model  due  to  Srinivasan  (Ref.  1) .  There  are  six  options  regarding  the  grid  generation 
schemes  employed  and  the  code  is  fully  compatible  with  the  QUICK  geometry  modeling  system.  The 
purpose  of  this  document  is  to  describe  the  code  in  a  manner  that  will  familiarize  the  user  with  the 
code’s  options  and  structure  and  enable  relatively  painless  usage  of  the  code. 

2.0  GENERAL  THEORY 


2.1  The  governing  equations 

The  governing  equations  solved  by  EULER3DS  are  integral  equations  representing  conservation  of 
mass,  momentum,  and  energy.  These  equations  are  now  described.  To  begin  with,  we  introduce  the 
essential  nomenclature.  The  dyadic  E  is  comprised  of  the  three  Cartesian  flux  vectors  and  is  given 
as 


F=ie+jf+kg 

where  the  vectors  e,f  and  g  are  given  by 


QU 

QV 

~  QW 

p  +  Qu} 

QUV 

QUW 

QUV 

f= 

p  +  QV^ 

8  = 

QVW 

QUW 

QVW 

p  +  QW 

_(e  +  p)u_ 

_(e+p)v_ 

_(e  +  p)w_ 

The  unit  vectors  i,y,  and^  are  the  basis  directions  in  the  Cartesian  x,y,  and  z  coordinate  system  and 
u,  V,  and  w  are  the  corresponding  velocity  components,  p  is  the  static  pressure  and  p  is  the  fluid  den¬ 
sity.  e  is  the  total  energy  per  unit  volume  given  by  the  equation 
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e  =  Q(ei  +  .5(m^  +  +  w^)) 


where  is  the  specific  internal  energy.  The  governing  integral  equation  is  then  given  as 


J(V*£)  dV=0 

V 

This  equation  is  transformed  into  a  boundary  integral  equation  using  Green's  theorem.  The  result 
is 


^n*F_dA  =  0  (2.1) 

dV 


In  this  equation,  ft  is  the  outward  unit  normal  vector  to  the  boundary  of  the  volume  V.  This  bound¬ 
ary  integral  equation  is  discretized  and  solved  on  each  of  the  finite  volume  cells.  This  discretization 
process  is  discussed  in  the  next  section. 

It  should  be  noted  that  in  the  actual  solution  algorithm,  only  four  equations  are  solved  since  for 
the  flows  of  interest,  the  total  enthalpy  is  constant.  However,  this  fact  is  not  used  in  the  development 
of  the  solution  algorithm,  thus  making  the  algorithm  more  general. 

2.2  The  Space  Marching  Philosophy 

For  this  discussion,  consider  the  general  six  sided  region  shown  in  Figure  1. 


I 


o  (rj=  constant) 
constant) 

^  constant) 
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Equation  (2.1)  must  be  applied  to  a  general  region  of  this  type.  We  first  recognize  that  the  quantity  MA 
is  locally  given  by 


MA  =  (r|  X  r,j)d^dT] 

on 

MA  =  {Ftj  X  r^)d$dt] 

on 

r 

X 

II 

on 

r 

ndA  =  (fj  X  f^)dt]d^ 

on 

r 

MA  =  (f^  X  r^dtd^ 

on 

MA  =  (f|  X  f^)dtd$ 

on 

v 

where  r  is  the  local  position  vector.  At  this  point,  we  choose  ?  to  be  the  marching  direction.  This 
means  that  S  =  constant  surfaces  are  solution  surfaces  and  we  must  solve  for  the  solution  surface  in 
terms  of  everything  else.  This  requires  that  the  component  of  Mach  number  in  the  C  direction  be  super¬ 
sonic.  We  then  restrict  these  ?  =  constant  surfaces  to  be  planar  and  proceed  to  approximate  the  inte¬ 
grals  over  them.  With  reference  to  Figure  2,  a  second  order  evaluation  on  the  surface  is  given  by 


/  X  rr,)*F  d^drj  ==  X  x  r^(D]  •  Fc  d^drj  (2.2) 


where  £c  is  the  Cartesian  flux  dyadic  at  the  quadrilateral  centroid.  The  actual  centroid  location  is 


Figure  2.  A  general  quadrilateral  region  on  =  constant 
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determined  as  follows.  For  the  general  planar  area,  S,  the  centroid  is 


rc  = 


When  this  formula  is  applied  to  the  quadrilateral  in  Fig.  2,  we  write 


or 


or 


rc  = 


jm 

f  rdSA  + 1 

1  dS 

.jd5 

where  Sa*  Sb  represent  the  areas  of  triangles  A  &  B  and  represent  the  centroids  of  these  tri¬ 

angles.  Now  the  triangle  centroids  are  given  simply  as 

=  j(^l  +  r2  +  f4) 

and  rB  =  jir2  +  r3  +  f4) 

where  ^1,2, 3, 4  represent  the  position  vectors  to  the  corners  of  the  quadrilateral  region.  The  right 
hand  side  of  Eq.  (2.2)  is  simply  expressed  as 


E  =  =  M 


where  C  is  a  unit  vector  perpendicular  to  ?  -  constant  surfaces  and  A  is  the  area  of  the  quadrilateral. 
The  vector  ?  is  given  by 
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Next,  consider  a  I  =  constant  surface.  Figure  3  illustrates  this  situation. 


This  surface  is,  in  general, 


non-planar.  Consequently,  its  contribution  to  the  integral  equation  (2.1)  is  slightly  different  than  that  for 
the  S  =  constant  surface.  We  seek  to  approximate 

\  (rti'X  t\)*Fdt]d^ 

r 

on  this  non-planar  quadrilateral.  We  define  the  quadrilateral  centroid  in  just  the  same  manner  as  previ¬ 
ously  done  for  the  C  =  constant  surface.  Thus,  the  integral  is  approximated  by 


F  =  ^•Fc 


where 


l  =  [|-[r,(n  X  +  +  X  h{rr)^dridt 

Note  that  the  bracketed  terms  in  this  expression  represent  the  MA  for  the  triangles  A  &  B  of  Figure  3. 
This  equation  is  simplified  to 
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^  =  X  ri;{r}^)l  +  [r^i^~)  x  r^iv~)'\  +  lr^iv^)  x  r^iv'Wtjdt 


Finally,  for  an  V  =  constant  surface,  we  define  the  appropriate  integral 

j  (ff  X  f|)  •  £  dl; 

t 


to  be 


G  =  »;•£ 

where 


X  f5(n]+ w)  X  x 


The  approximation  to  the  integral  equation,  Eq.  (2.1),  is  then  given  for  the  finite  volume  shown  in 
Figure  4  as 

=  M-)  -  [Rr)  -  Kt)]  -  m)  -  ^v)]  (2-3) 


Figure  4.  Generic  finite  volume  approximation 

from  which  it  is  clear  that  the  value  of  E  at  the  next  £  =  constant  surface  is  computed  at  each  integration 
step.  This  vector  E  must  then  be  decoded  to  yield  the  primitive  variable  information  for  the  next  step. 
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2.3  Decoding 

The  general  decoding  procedure  applies  to  the  vector  S  *  ^  •  We  consider  that  we  have  taken 

an  integration  step  and  that  A?*)  is  known.  We  define  the  basis  vector  ?  to  be  5  =  (Ci  .?2  .?3)  and 
further  define  the  contravariant  velocity,  W,  to  be 

where  q  is  the  fluid  velocity  vector,  q  =  (u,v,w).  We  then  write 


Ex 

E2 

£3 

£4 

Es 


=  qW 

=  p^x  +  QuW  =  p^i  +  uEx 
=  P^2  +  pvW  =  p^2  +  vEx 
=  P^3  +  =  pti  +  wEx 


-H.E, 


(2.4) 

(2.5) 

(2.6) 

(2.7) 


where  h  is  the  gas  enthalpy  and  Ht  is  the  constant  total  enthalpy.  We  define  the  ratio  of  enthalpy  to 
internal  energy  to  be 


Thus, 


yp _ ypW 

(f  - 1)0  "  (f  -  l)£i 


(2.8) 


Equations  (2.4)  through  (2.6)  are  solved  for  velocity  components  yielding 

M  =  (£2  -p  ^l)/Ex 

V  =  (E3  -p^2)/Ex  (2.9) 

W  =  (£4  -  p  ^3)/ £1 

Substitution  of  these  expressions  into  Eq.  (2.8)  yields 

h  =  ,_:^^-[Ci£2  +C2S3  +e3£4  -K?1  +  ?2  +  ?3)] 
kY  - 
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Substitution  of  Eqs.  (2.9)  and  (2.10)  into  Eq.  (2.7)  yields  a  quadratic  equation  for  pressure,  p. 
{2(7-^)  +^3£4]|p 

+  1^?  Ht  -y(ii  +El  +4  )j  =  o 

This  closed  form  solution  for  pressure,  P,  is  possible  only  in  the  case  of  7  =  constant  for  a  perfect 
gas.  From  this  value  of  P,  the  density,  p  ,  is  found  from  the  equation 


_ 

E2  Cl  +  £3  C2  +  E4  C3  -p(Ci  +  C2  +  C3) 


(2.11) 


which  is  obtianed  from  a  linear  combination  of  Eqs.  (2.4),  (2.5),  and  (2.6).  The  velocity  compo¬ 
nents,  u,  V,  and  w  are  determined  from  Eq.  (2.9). 

For  a  real  gas,  an  iterative  procedure  is  employed.  The  first  iteration  uses  7  to  establish  the 
initial  guess.  The  function  Q(p)  is  established  to  be 


Q(p)  =  h  + 


(2.12) 


A  Newton  iteration  yields 


pm+l  ^  pm 


Q'(pn 


where 


!2'(P) 


dp 


du  dv  dw 

u — +  V— +  w — 
dp  dp  dp 


(2.13) 


The  quantities 


du  dv  dw 
dp  *  dp  '  dp 


are  obtained  by  differentiating  Eqs.  (2.9). 


dp 


is  obtained  by 


differentiating  Eq.  (2.10).  This  latter  step  involves  the  constant  7  .  This  influences  the  convergence 
path  of  the  Newton  iteration  but  not  the  result.  The  general  iterative  procedure  then  is  to  initialize  p, 
p,  u,  V,  w,  using  the  constant  7  .  h  is  determined  by  the  real  gas  model  as  h  =  h(p,p).  Equation 
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(2.12)  gives  a  value  of  Q  and  Eq.  (2.13)  gives  Q’.  The  update  is  performed  generating  a  new  value 
of  p  from  which  Eqs.  (2.11)  and  (2.9)  give  p,  u,  v,  w.  This  cycle  continues  until  p  ceases  to  vary. 

2.4  Shock  Fitting 

EULER3DS  is  applicable  to  problems  in  which  the  Mach  number  in  the  I  -  direction  is  supersonic. 
Usually,  external  flow  problems  of  this  type  have  an  outer  shock  boundary.  As  a  result,  the  code  includes 
a  shock  fitting  algorithm  which  is  now  described. 

The  shock  is  assumed  to  be  a  g  =  constant  surface  which  separates  a  uniform  freestream  flow  from  its 
"shocked”  state.  The  g  =  constant  surface  is  taken  to  lie  just  downstream  of  the  shock  so  points  on  this 
surface  are  discontinuously  related  to  the  freestream  through  the  Rankine-Hugoniot  jump  conditions. 
Consider  the  shock  surface  in  Fig.  5. 


Figure  5.  Typical  g  =  constant  shock  surface 

Assume  conditions  ahead  of  the  shock  (eg.  :  pi,  qu  5i)  are  known  and  that  pressure,  Pi  ,  behind 
the  shock  surface  is  known.  The  conservation  equations  for  the  fluid  states  across  the  shock  are  given  as 


U\j>  “  U2j' 


where  subscripts  N,T  denote  normal  and  tangential,  respectively.  Introducing  yi  ,  yx  to  be 
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yi  =  —  ,  72  = 


hi 


the  following  equation  is  obtained  by  combination  of  the  above  conservation  laws. 


ul^  = 


(yi-i)(Pi-P2)  1 

'(y2  +  i)/i2  +  (y2-i)pi'l 

2qi  I 

AY2  -  l)Pl  -  (fl  -  1)^2] 

(2.14) 


Combination  of  continuity  and  normal  momentum  also  yields  an  equation  for  density,  Qi  . 


Q2  = 


Pi-Pz  +  Qiul^ 


(2.15) 


Thus,  for  a  perfect  gas  fi  =  72  =  y  and  Eq.  (2.14)  yields  when  P2  is  known.  For  a  real  gas,  the 
value  of  fz  is  unknown  inEq.  (2.14).  Consequently,  an  initial  guess  is  made  for  72  and  Eq.  (2.14)  yields 
a  value  of  Ui^  .  Equation  (2.15)  yields  Qi  .  The  combination  of  Pi  ,  Qi  yields  ht  ,  from  which  a 
new  Yi  is  determined.  A  repeat  cycle  yields  a  new  value  of  Qi  ,  etc.  Ths  is  repeated  until  Qi  converges. 

Once  this  portion  is  completed,  values  of  piy  Qi,  72,  /i2»  ^(2*  are  known.  The  value  of  ui^  is 
now  used  to  determine  how  to  propagate  the  shock  point.  In  general,  for  the  shock  shown  in  Fig.  5, 

rtj  X  ft 
n  =  — ^ ^ 

I  X  I'e  I 


and 


«ljv  =  9l  •  « 


(2.16) 


Since  generally  at  any  marching  station,  r,,  is  known  and  is  known,  the  vector  is  the  quantity  of 
concern  in  this  equation.  Equation  (2.16)  represents  a  single  scalar  equation  with  the  three  un¬ 
knowns  .  However,  zg  controls  the  way  in  which  the  physical  axial  coordinate,  z,  is  related  to 

the  computational  marching  coordinate,  ?  ,  and  as  such,  is  specified.  =  1  in  EULER3DS  and  conse¬ 
quently  are  the  two  unknowns.  An  additional  relationship  is  needed  in  order  to  solve  for  these 

quantities.  This  additional  relationship  controls  the  manner  in  which  a  shock  point  propagates  relative  to 
the  existing  shock  surface.  We  have  freedom  in  choosing  exactly  how  to  propagate  a  shock  point.  There 
are  three  options  coded  in  EULER3DS.  They  are: 


Option  1:  motion  along  existing  t|  =  constant  lines 

Option  2:  along  a  projection  of  ^  onto  a  ^  =  constnrit  'surface 

Option  3:  along  a  cylindrical  ray 
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These  options  are  now  described. 

Option  1:  In  this  case,  that  part  of  excluding  the  part,  say  ,  must  be  aligned  with  . 


Consequently,  we  may  write  .  Substitution  of  this  into  Eq.  (2.16)  yields 

one  equation  with  the  single  unknown,  s.  Solution  for  s  then  yields  the  shock  slope  vec¬ 
tor,  Tg  .  The  shock  position  may  be  updated  by  integrating  this  vector. 

Option  2:  In  this  case  is  expressed  as 


ft;  =  s{rr,  X  k) 


Again,  s  is  the  only  unknown  from  which  the  shock  slope  is  determined. 
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Option  3:  In  this  case,  h  is  given  by  =  ■s’(cos<^  /  +  sin0y)  where  <p  is  the  meridional  angle  of 

the  shock  point  in  a  cylindrical  system  defined  by  x  =  rcos(j)  ,  y  =  /"sin^  .  Again,  s  is 
the  only  unknown  from  which  the  shock  slope,  i  is  determined. 

The  shock  fitting  algorithm  thus  assumes  a  value  of  Pi  behind  the  shock.  The  remainder  of  the 
flow  variables  and  the  vector  are  determined  in  this  routine. 

It  should  be  noted  that  in  the  event  that  pi  <  p\  at  any  given  point  during  the  calculation,  the 

code  sets  Pi  -  p\  and  constructs  in  the  normal  way.  The  effect  of  this  is  that  if  a  shock  degener¬ 
ates  to  a  Mach  surface,  the  solution  is  maintained  correctly  and  the  point  is  propagated  along  the 
Mach  surface  correctly. 


3.0  DISCRETIZATION 

The  equatons  and  ideas  discussed  in  Secton  2.0  are  implemented  on  a  discretized  domain.  The 
purpose  of  this  section  is  to  quantify  this  domain  and  quantify  the  numerical  implementaion  of  the 
governing  equations  on  this  domain. 

3,1  Geometric  Related 

First  of  all,  the  code  assumes  pitch  plane  symmetry.  As  mentioned,  the  solution  to  the  governing 
flow  equations  is  obtained  by  marching  in  the  g  (i.e.,  z)  direction.  Thus,  a  discretized  grid  or  mesh  is 
required  at  each  g  marching  station.  This  mesh  has  as  its  four  boundaries  the  following 


^min  ^  vehicle  surface 
^max  ^  shock  suffacc 
timin  wind  symmetry  boundary 
Vmax^  symmetry  boundary 
The  actual  mesh  layout  is  now  described. 

3.1.1  Mesh  Layout 

The  'corner  grid'  is  defined  to  be  that  net  of  points  located  at  the  corners  of  the  finite  volume 
cells.  The  'cell  center  grid’  is  that  net  of  points  occupying  the  cell  centroid  locations.  The  coor¬ 
dinates  of  both  of  these  grids  are  known.  However,  the  flow  solution  is  generally  known  only  on  the 
'cell  center  grid'  except  at  boundaries. 

Since  the  integration  is  performed  in  the  ^  direction,  a  two-dimensional  grid  is  required  at  each 
^-station.  Figure  6  illustrates  the  typical  grid  layout  for  this  code. 
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Figure  6.  Mesh  Layout 


»  J 

corner  grid 
cell  center  grid 


Notice  that  the  centroid  grid  points  exist  on  the  symmetry  boundaries,  but  the  corner  grid  is  re¬ 
flected  across  these  boundaries.  Notice  also  that  cell  center  grid  points  are  defined  on  the  body  and 
the  shock.  This  is  to  provide  flow  data  on  these  surfaces  instead  of  one  half  cell  away.  These  points 
are  actually  not  cell  centers  but  rather  segment  centers.  The  coordinates  of  these  points  are  simply 
the  average  of  their  two  neighbors.  Figure  6  indicates  that  the  index  j  is  associated  with  the Tj  coor¬ 
dinate  and  k  is  asssociated  with  the  g  coordinate.  The  index  n  is  used  for  the  £  coordinate.  The 
’corner  grid’  consists  of  JX  by  KX  points  including  the  reflected  symmetry  boundary  points.  This 
means  that  the  'centroid  grid’  contains  JXl  by  KXPl  points  where  JXl  =  JX  -1  and  KXPi  =  KX  +  1 
There  are  JXl  by  KXl  finite  volume  cells  (KXl  =  KX  -  1).  The  variables  NTOT,  NCATOT,  and 
NCTOT  are  used  in  the  code.  They  are  defined  as 


NTOT  =  JX  *  KX 
NCATOT  =  JXl  ♦  KXl 
NCTOT  =  JXl  ♦  KXPl 


Thus,  there  are  NTOT  'corner  grid’  points,  NCATOT  cells,  and  NCTOT  ’centroid  grid’  points.  The 
flow  solution  is  computed  and  stored  at  the  NCTOT  ’centroid  grid’  points.  The  geometric  information 
is  stored  as  follows  : 
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XN,  YN  NTOT  *  corner  grid’  points  at  marching  station,  n. 

X  ,  Y  NTOT  'corner  grid'  points  at  marching  station,  n+1. 


XCN,  YCN  NCTOT  'centroid  grid*  points  at  n. 

XC  ,  YC  NCTOT  'centroid  grid'  points  at  n+1. 


VCELLN  NCATOT  cell  areas  at  n. 
VCELL  NCATOT  cell  areas  at  n+1. 


3.1.2  Grid  Generation 

The  corner  grid  is  determined  by  a  general  grid  generation  algorithm  with  several  options.  These 
options  include  the  scheme  of  Winslow  (Ref.  2),  Middlecoff’s  (Ref.  3)  application  of  Thompson’s 
(Ref.  4)  scheme,  Noack's  algebraic  scheme  (Ref.  5),  a  simple  algebraic  scheme,  and  Sorenson’s  (Ref. 
6)  application  of  Thompson's  scheme.  All  of  these  grid  generation  options  require  a  specification  of 
both  body  surface  points  and  shock  points.  The  body  surface  points  are  determined  by  the  QUICK 
geometry  modeling  system  (Refs.  1-9)  and  the  shock  points  are  determined  by  integration  of  shock 
slopes  determined  in  the  shock  fitting  scheme  discussed  in  section  2.4.  Once  the  corner  grid  is 
known,  the  cell  center  grid  is  determined  by  constructing  the  coordinates  of  all  of  the  quadrilateral 
centroids  as  described  in  section  2.2. 

3.1.3  Body  Surface  Description 

EULER3DS  is  designed  with  only  one  subroutine  which  interfaces  the  code  to  the  particular  ge¬ 
ometry  being  solved.  BPOINT  takes  as  input  the  value  of  z  (or  ^  )  and  returns  JX  values  of  x,y 
which  represent  the  body  cross  section  coordinates  at  this  z  station.  The  meridional  placement  of 
these  grid  points  on  a  given  body  cross  section  is  based  on  equal  arc  length  in  the  t|  direction. 
BPOINT  is  interfaced  to  the  QUICK  geometry  modeling  system.  For  the  given  z  value,  BPOINT 
breaks  the  cross  section  contour  into  NPMAX  points  (NPMAX  presently  equals  1000)  spaced  at 
equal  meridional  angle  increments.  The  coordinates  of  these  sample  points  in  this  fine  database  are 
obtained  by  calling  the  QUICK  subroutine  CSGEOM.  BPOINT  then  interpolates  the  JX  grid  points 
spaced  according  to  equal  arc  length  onto  this  fine  database  thus  providing  the  body  surface  'corner 
grid'  points. 

3.2  Integration  Algorithm 

This  section  expands  on  the  details  of  the  novel  integration  algorithm  used  in  EULER3DS.  We 
begin  with  the  generic  integration  scheme  applied  to  Equation  (2.3)  (repeated  here  for  convenience). 


m = m  -  m  -  m]  -  [W)  -  ^n)]  o-  d 


3.2.1  The  Generic  Integration  Scheme 

Consider  Figure  7  as  a  reference.  Equation  (3.1)  is  written  in  discretized  form 
as 


pn+1 


^M/2,k 


An+lf2 


f2 

^;,/c+i/2 


) 


(3.2) 
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Figure  7.  Discretized  generic  finite  volume  cell 

where  1  <  j  <  jxl  and  1  <  k  <  kxl.  This  equation  thus  produces  a  solution  at  all  cell  center 
points.  The  body  and  shock  points  are  treated  separately.  This  scheme  obviously  cannot  be  imple¬ 
mented  without  some  decision  about  how  to  evaluate  the  F  ,  G  quantities  at  the  cell  face  centers. 

Consider  first  the  $  =  constant  cell  face  on  which  we  wish  to  evaluate  F^jUjl k  ^)’ 


Figure  8.  Typical  ^  =  constant  cell  face 

From  this  sketch,  it  is  clear  that  the  four  closest  centroid  points  at  which  the  flow  data  are  known  are 

those  points  at  (j,k,n),  (j,k,n+l),  (j,k+l,n),  and  (j,k+l,n+l).  The  value  of  Fj^il2,k  obtained 

from  data  at  these  four  points.  Since  the  data  at  the  two  (n+1)  points  is  not  really  known  apriori, 
estimates  which  come  from  a  predictor  step  are  used.  This  makes  the  overall  scheme  a 
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predictor-corrector  scheme.  The  rationale  for  choosing  which  of  the  four  points  to  use  in  comput¬ 
ing  based  upon  the  signal  propagation  mechanisms  within  the  fluid.  This  subject  is  now 

discussed. 


3.2.2  Signal  Splitting 


The  vector  F  is  a  linear  combination  of  e,  f,  and  g  and  is  written  as  or 


F=|ie  +  |2/'+^3g 


The  vectors  e,  f,  and  g  may  be  written  as 

e  =  AU  ,  f=:BU  ,  g=CU 

where  U  =  (  p  ,  pu  ,  pv  ,  pw  ,  e  )  and  A,  B,  and  C  are  rotation  matrices  which  coincide  with  the 
de  df  dg 

Jacobian  matrices  dU  '  dU  '  dU  ,  respectively,  for  the  case  of  a  perfect  gas.  For  a  real 

de 

gas,  A  ^  ^  etc.,  but  e  =  AU  still  holds  for  the  appropriate  choice  of  A.  This  choice  is  simply 

h 

the  Jacobian  A  for  a  perfect  gas  with  y  replaced  by  ”7 .  Thus,  if  we  define  a  matrix  B  as 

5*  =  |,A  +  |25  +  |3C 

then  we  may  write 


F=B*U 

In  a  like  manner,  we  may  obtain 

E  =  A*U 


(3.3) 


(3.4) 


where  A*  =  5iA  +  $2^  +  SsC  and  C  =  (?i  Combination  of  Eqs.  (3.3)  and  (3.4)  yields 

F^B\Ay^E 

Since  E  is  the  dependent  variable  vector  and  the  flux,  F  ,  is  expressed  as  a  rotation  of  E  through 
the  diagonalizable  matrix  B* A*  F  may  be  split  into  components  according  to  the  system  character¬ 
istics  (eigenvalues  of  B* A*  ^).  In  fact,  B* A*  ^  may  be  written  as 


=  TKT^ 

where  A  is  a  diagonal  matrix  containing  the  eigenvalues,  A/  ,  of  B* A*  T  is  a  matrix  whose  columns 
contain  the  right  eigenvectors,  rt  ,  of  B* A*  ^ ,  and  7^^  is  a  matrix  whose  rows  contain  the  left  eigen- 

*  *~i 

vectors,  li  ,  of  B  A 

The  eigenvalues  of  B*A*  ^are  given  by 
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2,3  =  ' 


-[J'V  -  c2(e*|)  ±  c/MD 

/MZ) = c2(M)2  -  2-u -V  (M) + TJ  2(1-.  f) + -V  .  r)  -  c2(r*  r)(f  *  f) 

where  TJ  =  (^  •  ^  ,  "V  =  (§  •  ^  c  is  the  sonic  velocity.  The  left  eigenvectors,  h  ,  are  ob¬ 
tained  by  solving 

iB\Ay^-Xilfli  =  0  (3.5) 

This  is  a  formidable  task.  A  simpler  system  to  solve  is  obtained  by  recognizing  that 

A  =  MAM-^ 

de 

where  M  represents  a  similarity  transformation  and  A  is  a  Jacobian  matrix  where  U  is  the 

primitive  variable  vector.  Likewise, 


B  =  MBM^^  and 

C  =  MCM-^ 

from  which 

A*  =  MA*M-^ 

B*  =  MB*M-^ 

where 

A  =^iA4‘^2^+C3C  and 

r  =  ^iA  +  |2B  +  l3C 

Thus  we  write 

from  which  Eq.  (3.5)  becomes 

(5*-Aiiy/,  =  0 

where  if  ~  M^li  .The  vectors  k  are  more  easily  determined  and  are  found  to  be 


i/  =  («iic^  .  0  ,  aisc^  ,  -ai2C^  ,  -aj^) 
y  =  (-aiiai2C^  .  «i3C^  .  0  .  ,  «iiai2) 

y  =  («n«i3C^  ,  «i2C^  .  -«nc^  ,  0  ,  -anfliia) 

1/  =  (0  .  -Qc^a^i ,  -Qc^a^2  ,  -QC^a43 ,  "V  -I4TJ  ) 

=  (0  ,  -QC^asi ,  -Qc^a52  ,  -gc^ass ,  W  ) 
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The  matrix  M  is  well  known  to  be 


From  this  //  may  be  determined.  The  aj^  in  the  //  vectors  is  defined  to  be  Kk  • 

The  right  eigenvectors,  n  ,  are  obtained  by  solving 

(5*(A*)-'-AiI)r;  =  0 

Again  a  simpler  system  is  obtained  as  before  resulting  in 

from  which  r;  =  M  A*f/  .  The  are  found  to  be 


=  ,  0  ,  ai3 ,  -ai2  ,  0) 

h'^  =  —^(.-ai2  ,  aiittia ,  0  ,  -a?i  ,  0) 

aiiD 

^3^  =  — ^(«13  ,  «iiai2  ,  -ail  .0,0) 

anD 


Q  +  ai2  +  0:43) 
“V 


-0:41  ,  -a42  ,  -043  .  Q  ("V  -  X^  ) 


Q  {a\i  +  a52  +  Qsa) 

“v  -;i5'u 


-asi ,  -<252 


e  cv  -  W  ) 


where  D-~U  +  aj2  +  ^13) 

and  Dls  =  2g  c^ETJ  (^•F-X,J»  f)  -“V  (i>  C-X^.sC- 01 

From  this,  the  r;  are  determined. 

The  diagonalization  of  the  matrix  5*A*  ^allows  us  to  write 


22 


AEDC-TR-87-36 


5 

F=Fi  +  F2  +  F3  +  F4  +  F5=  ^  F 

i  =  1 

where  Ff  =  TAiT^^E  and  At  is  zero  everywhere  except  in  position  (i,i)  where  the  element  Xi  exists. 
With  some  manipulation,  it  is  easily  shown  that  this  simplifies  to 

Fi^XiriljE  (3.6) 


which  is  interesting  since  it  shows  that  the  flux  vector  component,  Fi  ,  has  the  direction  of  the 
right  eigenvector.  Since  A/  ^  ,  Eq.  (3.6)  is  also  written  as 

Fi  =  n-  If  F 

5 

The  significance  of  writing  F  =  ^  A/  A*/  /f  £  is  that  each  characteristic  piece  of  F  is  identified  and 

i  =  1 

the  evaluation  of  these  pieces  can  take  into  account  the  signal  propagation  directions  for  these  flux 
contributions. 

We  began  with  the  problem  of  determining  F  on  the  ^  =  constant  cell  face.  We  now 
have  the  problem  of  determining  five  values  of  Ff  on  this  face.  Two  candidate  evaluations  of  the  Fi 
are  introduced  in  the  following  two  sections. 

3.2.3  The  M-Flux 

The  M-flux  is  styled  after  a  finite  volume  version  of  the  MacCormack  scheme  as  applied  to  each 

5 

of  the  Ff  .  We  will  denote  this  flux  as  £m  =  ^  .  With  reference  to  Figure  8,  we  define  A/  to 

i  =  1 

be  Xi  calculated  using  the  ^  ?  values  corresponding  to  the  cell  face  under  consideration  and  the 

primitive  variables  on  that  face  which  are  obtained  by  averaging  the  nearby  neighbors.  Thus 

=  +  I".  ?)] 

where  U  is  the  primitive  variable  vector.  Also, 


a/  = 


Note  that  A/  is  an  effective  eigenvalue  on  the  g  =  constant  cell  face.  Then  we  write 


(3.7) 


Thus,  if  the  eigenvalue  is  positive  (  A/  >  0  ) ,  the  information  for  Ff  on  the  ^  =  constant  cell 
face  comes  from  points  "upwind”  at  n  and  "downwind”  at  n+1  or  along  the  characteristic.  The 
sgnlA/  1  coefficient  is  necessary  to  cause  Fm^  =  0  in  the  event  that  A/  =  0  as  indicated  by  Eq.  (3.6). 
The  quantity  evaluated  at  n+1  results  from  a  predictor  step  in  which 
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(3.8) 


In  practice,  the  evaluation  of  Eq.  (3.7)  is  not  unique  since  Fi  =Aj  r,-  if  E  .  Several  methods  of 
evaluation  were  tried  and  one  clearly  worked  best.  The  portions  of  Eq.  (3.6)  which  involve  the  geo¬ 
metric  data,  I  and  ff  ,  are  unique  and  therefore  result  in  no  ambiquity.  However,  the  parts  which 
involve  the  flow  variables  may  be  constructed  in  a  number  of  ways.  If  we  seek  evaluation  of  Ff  on 
the  g  =  constant  face  under  investigation  then  ^  f  are  clearly  known  at  level  n  and  are  known  for 
level  n+1  once  the  n+1  grid  is  established.  If  we  discover  that  A/  for  this  face  is  positive,  for  ex¬ 
ample,  then  we  evaluate  Fj  as 


Fi  =  XiriUE 

where  A;  =  A/  and  r/  if  E  is  evaluated  using  | ,  C  at  n  or  n+1  as  appropriate,  and  the  primitive  vari¬ 
ables  from  point  (j,k)  for  the  level  n  value,  and  point  (j,k+l)  for  the  level  n+1  value. 

It  must  be  noted  that  the  M-flux  evaluation  of  the  Ff  is  perfectly  consistent  with  a  characteristic 
like  method  even  though  the  n+1  quantity  is  ’downwind’. 

3.2.4  The  U-Flux 

The  U-flux  evaluation  candidate  results  from  a  belief  that  all  information  must  come  from  ’up¬ 
wind’  defined  in  the  following  sense:  If  A”  >  0  at  a  point  then  all  information  for  the  integration  step 
to  n+1  must  come  from  the  negative  coordinate  increment  side  of  this  point.  This  attitude  is  ridicu¬ 
lous  unless  X  =  0  but  never-the-less,  the  result  of  this  thinking  is  an  evaluation  of  Fi/.  as 


k+l-Oi 


where  subscript  e  denotes  ’extrapolated  value’  and 


(3.9) 


k+ai-sgn(Xi) 


(3.10) 


As  in  the  case  of  the  M-flux,  practice  has  shown  that  a  better  formula  is  to  use  the  appropri¬ 
ate  ^  ,  C  vectors  for  the  cell  face  under  consideration  and  use  the  primitive  variables  from  the  ’up¬ 
wind’  point  indicated  but  use  A/  instead  of  A/  at  the  upwind  point.  So,  for  instance,  in  Eq. 

(3.10)  is  evaluated  as 


(.Fi)jMOi  =  d/)(r/  Ij  ^jMOi 


where  the  left  and  right  eigenvectors  and  E  use  primitive  variables  from  point  (j,k+<7f)  and  ^  ^ 

values  appropriate  for  the  cell  face  in  question.  As  in  the  case  of  the  M-flux,  the  quantities  required 
at  n+i  are  assumed  to  be  known  via  an  upwind  predictor  step  in  which  Fi  is  obtained  by  Eq.  (3.8). 
Either  of  the  M-flux  or  the  U-flux  method  will  produce  stable  numerical  results.  These  two  schemes, 
however,  exhibit  opposite  phase  error  characteristics.  It  is  for  this  reason  that  these  two  schemes  are 
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combined  into  a  one  parameter  family  of  schemes  in  an  attempt  to  produce  a  scheme  superior  to 
either  of  the  individual  schemes.  The  result  is  called  the  ALPHA-scheme  which  is  the  subject  of  the 
next  section.  Fromm  (Ref.  10)  tried  something  similar  to  this  for  a  simple  scalar  linear  equation. 

3.2.5  The  ALPHA-scheme 

To  begin  the  discussion  of  the  ALPHA-scheme,  we  simply  define  the  flux  on  the  typical  g  =  con¬ 
stant  cell  face  to  be  a  linear  combination  of  the  M-flux  and  the  U-flux.  That  is,  we  take 


Fi  =  a  Fmi  +  (1  -  Fu- 

When  a  =  1,  we  use  entirely  the  M-flux  and  a  =  0  gives  the  U-flux.  A  value  of  a  =  .5  weights  each 
part  equally  and  exhibits  superior  dispersive  properties  to  either  M  or  U  separately.  Experience  has 
shown  that  for  either  M  or  U  separately,  dispersion  grows  dramatically  worse  as  the  local  Courant 
number  becomes  small.  The  a-scheme  with  a  =.5,  however,  is  quite  insensitive  to  Courant  number. 

This  section  is  ended  with  a  summary  of  the  overall  integration  algorithm.  Figure  9  is  used  as  a 
reference. 
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From  Eq.  (3.2),  we  write 


predictor  ; 


=  ^j,k+l  ~  1  [(ft?.  ifc+1  ”  ,ic+i  C^/Oy.fe+i]  (3.11) 


/  =  1 


i=  1 


where,  for  example,  on  the  ^  =  constant  face, 


with 


and 


and 


imk^sgnmumjMai 


It  is  understood  that  quantities  such  as  WOy.fe  are  evaluated  using  the  primitive  variables  at  solution 
point  (j,k,n)  and  values  of  | ,  C  on  the  face  in  question  (eg.,  in  this  case  ,  ^  ).  This  is  be¬ 
cause  $  ,  Z  are  face  quantities  and  have  no  meaning  at  solution  points. 

Thus,  application  of  Eq.  (3.11)  gives  as  a  predicted  solution.  The  decoding  routine  dis¬ 
cussed  in  section  (2.3)  will  provide  the  primitive  variables  at  this  solution  point  provided  a  value 

for  is  known.  This  requirement  essentially  translates  to  knowing  the  grid  at  n+1  which  is  ob¬ 
tained  by  the  grid  generation  scheme  once  the  body  and  shock  boundary  points  are  known.  The 
body  points  are  known  at  n+1  by  virtue  of  the  QUICK  geometry  modeling  system,  and  the  shock 
points  are  known  at  n+1  by  integrating  the  level  shock  slope  vector,  ,  obtained  in  the  shock 
fitting  routine.  Thus,  decoding  is  done  and  all  information  at  n+T  is  available.  The  corrector  step  is 
then  taken  using 

corrector. : 


5  r 


where,  for  example 
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with  given  by  Eq.  (3.7)  and  Fy,  given  by  Eqs.  (3.9)  and  (3.10).  Note  that  the  direction  deci¬ 
sion  depends  on  Oi  .  This  must  be  based  on  level  A/  values  both  in  the  predictor  and  the  cor¬ 
rector.  However,  in  the  evaluation  of  (F/)^,  such  as  is  needed  in  Eq.  (3.7)  the  A/  used  there  is 
obtained  from  n+1  data. 

3.3  Boundary  Conditions 

The  integration  algorithm  described  in  the  previous  section  is  applicable  to  general  solution 
points.  There  are  four  boundaries  which  require  special  treatment.  The  two  symmetry  boundaries  are 
handled  with  a  simple  reflection  procedure.  The  body  boundary  condition  influences  the  solution  at 
the  solution  point  just  off  the  body  since  computation  of  this  point  requires  the  body  flux  information. 
This  flux  information  simply  enforces  the  surface  tangency  requirement.  The  solution  at  the  point 
directly  on  the  body  is  determined  with  a  half  cell  matching  procedure.  Following  each  predictor  and 
corrector  step,  the  surface  tangency  condition  is  enforced  explicitly  with  Abbett's  method  (Ref.  11). 

At  the  shock  surface,  the  point  at  k  =  kx  just  a  half  cell  inside  the  shock  requires  the  flux  at  the 
shock.  The  algorithm  upwinds  normally  based  on  the  eigenvalues  except  that  when  data  is  needed 
outside  the  shock  cell  face  for  the  flux  at  this  face,  the  values  on  the  shock  are  used  instead.  The 
shock  point  is  computed  with  a  half  cell  matching  procedure  identical  to  that  at  the  body.  This  pro¬ 
vides  a  complete  solution  at  the  shock  from  which  only  the  pressure  is  kept.  The  shock  fitting  routine 
is  then  used  to  enforce  the  Rankine-Hugoniot  jump  conditions. 

4.0  THE  STARTUP  PROBLEM 

EULER3DS  solves  the  inviscid  steady  flow  equations.  Consequently,  the  initial  solution  surface 
from  which  it  starts  marching  must  theoretically  be  a  correct  steady  state  solution.  This  can  be 
achieved  in  one  of  two  ways  generally.  If  the  nosetip  of  the  body  in  question  is  a  conical  nosetip, 
then  an  ” arbitrary”  initial  solution  may  be  marched  with  EULER3DS  down  this  hypothetical  conical 
body  (continuation  of  conical  nosetip)  until  a  steady  conical  solution  is  reached.  This  conical  solution 
may  then  be  scaled  back  to  the  nosetip  station  and  used  as  the  initial  solution  surface  for  the  actual 
body.  The  second  way  applies  to  bodies  with  non-conical  nosetips.  In  this  case  a  separate  computer 
code  is  required  to  solve  for  the  steady  state  solution  as  the  time  asymptotic  solution  to  the  unsteady 
equations.  The  first  of  these  two  ways  is  the  only  method  utilized  in  this  work. 

4.1  Initial  Solution 

Before  EULER3DS  can  begin  computing,  an  initial  solution  msut  be  provided.  EULER3DS 
reads  this  initial  solution  from  FORTRAN  unit  TRSTRT’.  The  sequence  of  read  statements  used  is 

READ  (IRSTRT,  END=4  444)  JX,  KX 

READ(IRSTRT)XREFM,YREFM,ZREFM,XFOR,YFOR,ZFOR,XMOM,YMOM,ZMOM 

NTOT=JX*KX 

NCTOT=  (JX- 1 )  *  (KX+ 1) 

RE  AD  (IRSTRT)  (P(L),RHO(L),U(L),V(L),W(L),L=l,NCTOT) 

,(X(L),Y(L),L=l,NTOT),Z 

The  information  required  by  these  read  statements  must  have  been  generated  and  written  out  in  this 
manner  by  the  user  prior  to  execution  of  EULER3DS. 

The  first  read  statement  expects  the  values  of  JX,  KX.  These  represent  the  number  of  points  in 
the  *  corner  grid’  around  the  body  and  between  body  and  shock  respectively. 

The  second  read  statement  reads  information  concerning  force  and  moment  loads.  The  vector 
(XREFM,  YREFM,  ZREFM)  represents  the  coordinates  of  the  point  about  which  the  moment  due  to 
pressure  is  computed.  The  vector  (XFOR,  YFOR,  ZFOR)  represents  the  accumulated  force  vector 
for  that  portion  of  the  body  upstream  of  the  current  station.  This  vector  is  structured  as 
(-NORMAL,  SIDE,  AXIAL).  The  vector  (XMOM,  YMOM,  ZMOM)  is  the  accumulated  moment 
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about  (XREFM,  YREFM,  ZREFM)  due  to  the  pressure  forces.  Note  that  the  units  on  these  quanti¬ 
ties  depend  upon  the  units  of  the  code's  (x,y,z,p)  values.  The  code’s  p  values  are  really  p/PINF  and 
it’s  x,y,z  units  are  arbitrary.  This  results  in  a  force  per  unit  reference  area  per  unit  pressure  and  a 
moment  per  unit  reference  area  squared  per  unit  pressure.  These  are,  in  effect,  force  and  moment 
coefficients.  Note  that  the  reference  area  is  the  square  of  the  physical  length  of  a  unit  x,y,  or  z  code 
element. 

The  third  read  statement  reads  the  values  of  pressure,  density,  and  velocity  components  at  each 
point  in  the  ’cell  center  grid’  (NCTOT  points)  and  the  grid  coordinates  of  the  ’corner  grid’  (NTOT 
points)  plus  the  value  of  z  at  this  station.  See  Figure  6  in  section  (3.1.1)  for  a  reference  on  the 
mesh  layout. 

The  user  is  responsible  for  providing  this  information  initially  to  EULER3DS  to  begin  computa¬ 
tion.  EULER3DS  produces  as  output  on  FORTRAN  unit  ’IRSTRT  +  1’  subsequent  solutions  in  this 
format  for  future  startups. 


5.0  INPUT  DATA 

In  addition  to  the  initial  solution  data  file  on  FORTRAN  unit  ’IRSTRT’,  EULER3DS  also  reads 
input  from  unit  5  and  unit  15.  The  data  on  unit  15  is  geometry  definition  data  which  is  created  by 
the  QUICK  geometry  modeling  system  described  in  Refs.  (7-9) .  A  description  of  this  package  is  be¬ 
yond  the  scope  of  this  report.  FORTRAN  unit  5  is  read  by  EULER3DS  to  provide  various  parame¬ 
ters  which  control  various  aspects  of  the  code.  These  parameters  are  described  in  the  next  subsec¬ 
tion. 

5.1  Input  Parameter  Description 

The  input  parameters  required  by  EULER3DS  are  cast  into  two  basic  categories.  These  are 
’general  input’  and  ’grid  related  input’.  The  ’general  input’  supplies  parameters  via  the  FORTRAN 
NAMELIST  called  ’NAMLST’  while  ’grid  related  input’  is  supplied  through  a  FORTRAN  NAMELIST 
called  ’GRDNAM’.  The  parameters  contained  in  ’NAMLST’  are  now  listed  with  their  defaults  indi¬ 
cated  in  parantheses  and  a  brief  description  of  each  follows. 


IRSTRT 

(10) 

NGAS 

(0)  ■ 

lORDER 

(1) 

ISHOCK 

(3) 

NSTEP 

(1) 

NOUT 

(1) 

lOUT 

(0) 

NGRD 

(100) 

ALFA 

(0.5) 

CN 

(0.95) 

FSMACH 

(no  default) 

PINF 

(1.0) 

RHOINF 

(1.0) 

RELAX 

(1.0) 

GAM 

(1.4) 

ALPHA 

(0.) 

NPLT 

(10000) 

NDSK 

(10000) 

IRSTRT:  FORTRAN  unit  number  for  initial  data  surface  (NOTE  :  output  solution  is  on  unit  irstrt+1) 
NGAS:  (=  0  for  perfect  gas)  or  (=  1  for  equilibrium  air  real  gas) 
lORDER:  (=  1  for  first  order)  or  (=  2  for  second  order) 

ISHOCK:  (=  1  for  shock  point  propagation  along  existing  T|  =  constant  lines) 

or  (=  2  for  propagation  along  projection  of  shock  normal  onto  ^  =  constant  surface) 
or  (=  3  for  propagation  along  cylindrical  rays) 

NSTEP:  number  of  marching  steps  this  run 

NOUT:  print  out  frequency  on  unit  6  (print  every  NOUT  steps) 
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lOUT:  [0=abbreviated  print  out  (body-shock  data  only)]  [l=full  flowfield  printout] 

NGRD:  max  number  of  iterations  in  grid  construction 

ALFA:  alpha  scheme  parameter  (0.=  pure  upwind,  l.=  modified  Mac.)  not  used  when  iorder=l 
CN:  Courant  number,  max  value=l  for  stability 
FSMACH:  freestream  Mach  no. (no  default) 

PINFSI:  freestream  pressure  in  SI  units  (needed  for  real  gas  option  only!  do  not  include  for  perfect 
gas) 

RINFSI:  freestream  density  in  SI  units(needed  for  real  gas  option  only!  do  not  include  for  perfect 
gas) 

RELAX:  relaxation  parameter  used  in  grid  line  relaxation  scheme 
GAM:  perfect  gas  specific  heat  ratio 
ALPHA:  angle  of  attack  in  degrees 

NPLT:  plot  output  frequency  for  genplot  (plot  every  NPLT  steps) 

NDSK:  restart  solution  output  frequency  (write  out  every  NDSK  steps) 

Next,  the  parameters  contained  in  ’GRDNAM'  are  listed  with  their  defaults  indicated  in  paren¬ 
theses  and  a  brief  description  of  each  follows. 

MGRD  (4) 

EPSSP  (no  default) 

EPSCP  (no  default) 

MGRD:  (l=elliptical  Laplace’s  equation  grid) 

or  (2=ellipticai  Poisson  equation  (p/q  by  Ref  3.)) 
or  (3=algebraic  grid  with  Laplace  smoothing  (Ref  5.)) 
or  (4=simple  algebraic  grid  (equal  space  from  body  to  shock)) 
or  (5=elliptical  Poisson  equation  (p/q  by  Ref  6.)) 

EPSSP:  Used  only  when  MGRD=3 

(0.  =  no  orthogonality  at  wall) 
or  (1.  =  max  orthogonality  at  wall) 

EPSCP:  Used  only  when  MGRD=3 

(0.  =  no  clustering  toward  wall) 
or  (1.  =  max  clustering  toward  wall) 

NOTE:  For  option  5,  the  p/q  control  functions  are  based  on  a  requirement  of  orthogonality  at 
the  body  and  first  cell  spacing  of  the  initial  grid  estimate  before  relaxation.  This  may  not  be  satisfac¬ 
tory  in  all  cases.  The  requirement  of  orthogonality  is  met  by  bisecting  the  angle  between  (j-l,j)  and 
(j,j+l)  segments.  This  procedure  gives  reasonable  results  even  at  a  pointed  tip. 

The  only  parameters  one  should  usually  need  to  prescribe  to  run  a  perfect  gas  solution  are: 

lORDER  (if  second  order  is  desired) 

NSTEP 

NOUT 

FSMACH 

ALPHA 


6.0  SUMMARY 

EULER3DS  is  a  robust  3-D  space  marching  code  for  external  inviscid  supersonic  flows  about 
arbitrary  geometries.  The  code  has  a  variety  of  options  for  various  circumstances.  The  use  of 
IORDER=l  is  recommended  for  situations  with  exceptionally  strong  shocks,  or  steep  gradients,  particu¬ 
larly  during  sudden  startups.  Once  this  sudden  behavior  has  passed,  the  user  may  switch  to 
IORDER=2. 

Caution  must  be  used  in  choosing  the  parameter  ISHOCK.  One  helpful  feature  of  EULER3DS 
is  that  points  on  the  body  contour  are  spaced  based  on  equal  arclength.  This  maintains  a  uniformity 
in  the  body  grid.  The  shock  point  spacing,  however,  is  an  ever  evolving  situation  depending  on  the 
shock  shape  and  its  evolution  in  the  £  direction,  and  on  the  way  the  shock  points  are  propagated. 

The  method  with  the  greatest  chance  of  success  for  an  arbitrary  problem  is  ISHOCK=3  (i.e,,  propaga¬ 
tion  along  cylindrical  rays  in  a  £  =  constant  plane.  When  this  occasionally  causes  shock  points  to  be- 
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come  spaced  improperly  over  some  ^  distance,  the  user  must  explore  the  other  options.  The  author 
has  developed  a  method  which  looks  very  promising  for  propagating  shock  points  based  on  the  shape 
of  the  body  contour.  This  option  is  not  present  in  EULER3DS  however. 

Regarding  the  grid  generation  schemes,  the  best  advice  is  to  use  the  simplest  scheme  which  will 
work  for  the  particular  problem  being  solved.  The  order  of  preference  would  be  MGRD=4,3, 1,2,5. 
]VIGRD=4,3,or  1  will  work  for  most  problems.  MGRD=2  and  5  represent  a  level  of  sophistication 
which  results  in  clustering  capability  at  the  expense  of  robustness  for  some  shapes. 
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